#!/usr/bin/env python
# coding: utf-8

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
import ipywidgets as widgets
np.random.seed(seed=43)
get_ipython().run_line_magic('matplotlib', 'inline')


# # Metodologia

# ## Geração de retornos e preços por Movimento Geométrico Browniano (caso base) ou sistema de mudança de regime.
# 
# **MGB**
# $$
# r_t = \mu + \sigma*\sqrt{dt}
# $$
# 
# **Mudança de Regime**: $\Delta \mu$ representa uma alteração repentina de preço e $I^t_{\{0,1\}}$ é uma variável binária que quando 1 indica que ouve a mudança no período $t$.
# $$
# r_t = \mu +\Delta \mu*I^t_{\{0,1\}} +\sigma*\sqrt{dt}
# $$
# 
# 
# $$
# P_{t+1} = P_t*e^{r_t}
# $$

def gbm(n_periods = 52, n_scenarios=10000, mu=0.00, sigma=0.1703, steps_per_period=1, s_0=191.89):
    dt = 1/steps_per_period
    n_steps = int(n_periods*steps_per_period) + 1

    rets_plus_1 = np.random.normal(loc=(1+mu)**dt, scale=(sigma*np.sqrt(dt)), size=(n_steps, n_scenarios))
    rets_plus_1[0] = 1
    return s_0*pd.DataFrame(rets_plus_1).cumprod(),  pd.DataFrame(rets_plus_1-1)

def regimen_change(n_periods = 52, n_scenarios=10000, mu=0.00, sigma=0.1703, steps_per_period=1, s_0=191.89, delta_mu=0.5):
    t_star = np.random.randint(1, n_periods-1, (1, n_scenarios))

    dt = 1/steps_per_period
    n_steps = int(n_periods*steps_per_period) + 1

    rets_plus_1 = np.random.normal(loc=(1+mu)**dt, scale=(sigma*np.sqrt(dt)), size=(n_steps, n_scenarios))
    rets_plus_1[t_star, np.array(range(n_scenarios))] = (1 + delta_mu)**dt
    
   
    rets_plus_1[0] = 1
    return s_0*pd.DataFrame(rets_plus_1).cumprod(),  pd.DataFrame(rets_plus_1-1)


# ## Cálculo de default
# 
# * A função utilizada para calcular a porcentagem de default é $d_\% = \sqrt{\mid{r_t}\mid}$
# 
# 
# ## Posições em aberto na bolsa
# Baseado na BBCE (Jan 2020): 
# * 1419 novos contratos por semana
# * 3.965.000 MWh negociados semanalmente
# 
# ## Posições assumidas pela clearing em caso de default
# $$
# x_t = \sqrt{\mid{r_t}\mid}*3.965.000*side
# $$ 
#  
# \begin{equation}
#     side=
#     \begin{cases}
#       -1, & \ r_t > 0 \\
#       1, & r_t < 0
#     \end{cases}
#   \end{equation} 
#   
# ## Posições em aberto na bolsa
# Baseado na BBCE (Jan 2020): 
# * 1419 novoa contratos por semana
# * 1419*2 = 2838 contratos abertos (considerando as duas pontas)
# * 3.965.000 MWh negociados semanalmente          
#               

def default_position (ret, open_positions = 3965000):
    return ret.applymap(lambda x: -np.sqrt(abs(x))*open_positions/100 if x>0 else np.sqrt(abs(x))*open_positions/100)


# ## Posição corrente da clearing 
# Baseado no índice de liquidez $L$ que representa quantas semanas a clearing fica presa com a posição.
# $$
# X_t = \sum_{k=1}^{L} x_{t-k}
# $$

def cumulative_clearing_position (default, L):
    return default.shift(1).rolling(L,min_periods=1).sum()


# ## Resultado financeiro da clearing
# Em cada período $t$:
# $$
# \pi_t = (P_t - P_{t-1})*X_t
# $$
# 
# Acumulado até o período $T$:
# $$
# \Pi_t = \sum_{t=0}^T \pi_t
# $$

def wealth_change (prices, X):
    return (prices - prices.shift(1))*X
    
def cumulative_wealth (wealth):
    return wealth.cumsum()


# ## Opção real do seguro (put, strike = 0)
# payoff $\alpha$ da opção no tempo $T$ no cenário $s$
# $$
# \alpha_s = max(-\Pi_T^s, 0)
# $$
# 
# payoff médio nos $S$ cenários:
# $$
# \alpha = \frac{1}{S} * \sum_{s=1}^S \alpha_s
# $$
# 
# valor da opção $p$ no tempo $t=0$ baseado na taxa de desconto livre de risco $r_f$
# $$
# p = \alpha*e^{-T*r_f}
# $$

def average_payoff(cum_pl):
    return cum_pl.iloc[-1].apply(lambda x: max(-x, 0 )).mean()

def present_value(value, rate, years = 1):
    return value*np.exp(-1*rate*years)
    

# # Resultados

# ## Variáveis
#  * drift de retornos $\mu$: baixo(-.05), base(0), alto(+.05) 
#  * Coeficiente de liquidez $L$: 1,2,3,4 semanas
# 
#  * volatilidade ($\sigma$) 17,03% por semana
#  * Tempo de simulação: 52 semanas

def plot_wealth(wealth, title, L):
    y_max=wealth.values.max()
    terminal_wealth = wealth.iloc[-1]
    
    tw_mean = terminal_wealth.mean()
    tw_median = terminal_wealth.median()
    failure_mask = np.less(terminal_wealth, 0)
    n_failures = failure_mask.sum()
    p_fail = n_failures/len(terminal_wealth)

    e_shortfall = np.dot(terminal_wealth, failure_mask)/n_failures if n_failures > 0 else 0.0

    # Plot!
    fig, (wealth_ax, hist_ax) = plt.subplots(nrows=1, ncols=2, sharey=True, gridspec_kw={'width_ratios':[3,2]}, figsize=(24, 9))
    fig.suptitle("Resultado Financeiro Clearing ("+ title +" - L = "+str(L)+")", fontsize=16)
    plt.subplots_adjust(wspace=0.0)
    
    wealth.plot(ax=wealth_ax, legend=False, alpha=0.3, color="black")
    wealth_ax.axhline(y=0, ls="--", color="black")

    wealth_ax.set_ylim(top=y_max)
    
    terminal_wealth.plot.hist(ax=hist_ax, bins=50, ec='w', fc='black', orientation='horizontal')
    hist_ax.axhline(y=0, ls="--", color="black")
    hist_ax.axhline(y=tw_mean, ls=":", color="black")
    hist_ax.axhline(y=tw_median, ls=":", color="black")
    #hist_ax.annotate(f"Mean: ${int(tw_mean)}", xy=(.7, .9),xycoords='axes fraction', fontsize=24)
    #hist_ax.annotate(f"Median: ${int(tw_median)}", xy=(.7, .85),xycoords='axes fraction', fontsize=24)
    
    wealth_ax.set_title('Resultado Acumulado')
    hist_ax.set_title('Distribuição Lucro Final')
    hist_ax.axhline(y=0, ls="--", color="black", linewidth=3)
    hist_ax.annotate(f"Option exercised: {n_failures} ({p_fail*100:2.2f}%)\nE(payoff): ${average_payoff(wealth):,.2f}", xy=(.2, .15), xycoords='axes fraction', fontsize=24)


price_up, returns_up = regimen_change(delta_mu=0.5)
price_up.to_excel("price_up.xlsx", index = True)
returns_up.to_excel("returns_up.xlsx", index = True)

price_base, returns_base = gbm(mu=0)
price_base.to_excel("price_base.xlsx", index = True)
returns_base.to_excel("returns_base.xlsx", index = True)

price_down, returns_down = regimen_change(delta_mu=-0.5)
price_down.to_excel("price_down.xlsx", index = True)
returns_down.to_excel("returns_down.xlsx", index = True)


# ## Comparação entre as posições médias mantidas pela clearing para diferentes liquidez (L = 1 ... L = 4)


def average_position(returns, L):
    x = default_position(returns)
    X = cumulative_clearing_position(x, L)
    return X.mean(axis=1)


avg_neg = {'L1': average_position(returns_down, L=1), 
           'L2': average_position(returns_down, L=2), 
           'L3': average_position(returns_down, L=3),
           'L4': average_position(returns_down, L=4)}

pd.DataFrame(avg_neg).to_excel("avg_neg.xlsx", index = True)

ax_down = pd.DataFrame(avg_neg).plot.line(color="black",style=["-","--","-.",":"])
ax_down.set(xlabel='semana', ylabel='posição (MWh)',
title='Posição Média da Clearing (Choque de Baixa(-50%))')

avg_base = {'L1': average_position(returns_base, L=1), 
           'L2': average_position(returns_base, L=2), 
           'L3': average_position(returns_base, L=3),
           'L4': average_position(returns_base, L=4)}

pd.DataFrame(avg_base).to_excel("avg_base.xlsx", index = True)

ax_down = pd.DataFrame(avg_base).plot.line(color="black",style=["-","--","-.",":"])
ax_down.set(xlabel='semana', ylabel='posição (MWh)',
title='Posição Média da Clearing (Sem Jumps)')

avg_pos = {'L1': average_position(returns_up, L=1), 
           'L2': average_position(returns_up, L=2), 
           'L3': average_position(returns_up, L=3),
           'L4': average_position(returns_up, L=4)}

pd.DataFrame(avg_pos).to_excel("avg_pos.xlsx", index = True)

ax_down = pd.DataFrame(avg_pos).plot.line(color="black",style=["-","--","-.",":"])
ax_down.set(xlabel='semana', ylabel='posição (MWh)',
title='Posição Média da Clearing (Choque de Alta (50%))')


f, (ax1, ax2, ax3) = plt.subplots(3, sharex=True, sharey=True)
ax1.plot(pd.DataFrame(avg_neg),'-',color='black')
ax2.plot(pd.DataFrame(avg_base),'-',color='black')
ax3.plot(pd.DataFrame(avg_pos),'-',color='black')

ax1.set_title('Posição Média da Clearing \n Choque de Baixa (-50%)')
ax2.set_title('Sem Jumps')
ax3.set_title('Choque de Alta (50%)')
ax3.set(xlabel='semana')
ax2.set(ylabel='posição (MWh)')
plt.subplots_adjust(hspace = 0.5) 


# ## Valor da opção para diferentes L e drifts


def option_value (prices, returns, L, title):
    x = default_position(returns)
    X = cumulative_clearing_position(x, L)
    pi = wealth_change(prices, X).dropna()
    PI = cumulative_wealth(pi).dropna()
    pd.DataFrame(PI).to_excel('cumulative'+ str(L) + str(title) + '.xlsx', index = True)
    plot_wealth(PI, title, L)
    return(PI.iloc[-1].mean())


DownL1 = option_value(price_down,returns_down, 1, 'Mudança de Regime Negativa (-50%)')
DownL2 = option_value(price_down,returns_down, 2, 'Mudança de Regime Negativa (-50%)')
DownL3 = option_value(price_down,returns_down, 3, 'Mudança de Regime Negativa (-50%)')
DownL4 = option_value(price_down,returns_down, 4, 'Mudança de Regime Negativa (-50%)')
BaseL1 = option_value(price_base,returns_base, 1, 'Caso Base')
BaseL2 = option_value(price_base,returns_base, 2, 'Caso Base')
BaseL3 = option_value(price_base,returns_base, 3, 'Caso Base')
BaseL4 = option_value(price_base,returns_base, 4, 'Caso Base')
UPL1 = option_value(price_up,returns_up, 1, 'Mudança de Regime Positiva (50%)')
UPL2 = option_value(price_up,returns_up, 2, 'Mudança de Regime Positiva (50%)'  )
UPL3 = option_value(price_up,returns_up, 3, 'Mudança de Regime Positiva (50%)')
UPL4 = option_value(price_up,returns_up, 4, 'Mudança de Regime Positiva (50%)')

print('Down L1 = ' + str(DownL1))
print('Down L2 = ' + str(DownL2))
print('Down L3 = ' + str(DownL3))
print('Down L4 = ' + str(DownL4))
print('Base L1 = ' + str(BaseL1))
print('Base L2 = ' + str(BaseL2))
print('Base L3 = ' + str(BaseL3))
print('Base L4 = ' + str(BaseL4))
print('UP L1 = ' + str(UPL1))
print('UP L2 = ' + str(UPL2))
print('UP L3 = ' + str(UPL3))
print('UP L4 = ' + str(UPL4))


